
t = 75;


Nx = params.Nx;
Ny = params.Ny;
Nt = Nx * Ny;

% u = deval(sol, tV);
% T  = reshape(u(0*Nt+1:1*Nt), Nx, Ny);
% ux = reshape(u(1*Nt+1:2*Nt), Nx, Ny);
% uy = reshape(u(2*Nt+1:3*Nt), Nx, Ny);
% P1 = reshape(u(3*Nt+1:4*Nt), Nx, Ny);
% P2 = reshape(u(4*Nt+1:5*Nt), Nx, Ny);
% 
% imagecsc(ux);

format compact
fprintf('RelTol: %g    AbsTol: %g\n', params.RelTol, params.AbsTol);


fprintf('**** t = %.3f ****\n', t);
ti = find(sol.x>t, 1, 'first');
fprintf('ti = %d, t = %.3f\n', ti, sol.x(ti));
timestep_ref = sol.x(ti) - sol.x(ti-1);
fprintf('ode45 timestep_ref =  %.6f\n', timestep_ref);


tord = 4;
d1k = d1ker(tord);
d2k = d2ker(tord);
abserr_dt = [];
relerr_dt = [];
abserr_dtdt = [];
relerr_dtdt = [];
pwrs = [(1:-0.25/2:-2)];

u = deval(sol, t);
ud1 = hydro_eom2(0, u, params);
ud2 = hydro_deom2(0, u, ud1, params);

for i = 2.^pwrs
  timestep = timestep_ref * i;
  fprintf('ode45 timestep =  %.6f\n', timestep);
  t_eval = t + timestep.*[-tord/2:tord/2];
  t_eval

  u_eval = deval(sol, t_eval);
  ud1_fd   = conv2(u_eval, d1k', 'valid')./timestep;
  ud2_fd = conv2(u_eval, d2k', 'valid')./timestep./timestep;
  [a r] = cmpder(u, ud1, ud1_fd, params);
  %[a r] = cmpder(u, ud2, ud2_fd, params);

  abserr_dtdt = [abserr_dtdt; a];
  relerr_dtdt = [relerr_dtdt; r];
end

plot(pwrs, log(relerr_dtdt))




%imageh(params, u, 1, 'T')




return;
% % investigate error of finite differencing of udt in order to determine optimal time step
% fail_udt = 1;
% while (fail_udt > 0)
%   ts = [-tord:1:tord] .* timestep + tOut(ti)
%   % evaluate ``true'' udt
%   [u, udt] = deval(sol, ts);
% 
%   % compute finite difference udt
%   udt_fd = conv2(u, d1k', 'valid')./timestep;
%   err_udt = abs(udt - udt_fd);
%   fail_udt = nnz(err_udt>params.AbsTol) & (err_udt>params.RelTol .* abs(udt));
% 
%   if fail_udt > 0
%     time_step = 0.9 .* time_step;
%   end
% 
% udtdt2 = conv2(udt, d1k', 'valid')./timestep;
% 
% udt = udt(:, tord/2+1:end-tord/2);
% udtdt = udtdt(:, tord/2+1:end-tord/2);
% 
% Tdt  = reshape(udtdt(1:params.Nx*params.Nx,:), params.Nx, params.Nx, size(udt,2));
% Tdt2 = reshape(udtdt2(1:params.Nx*params.Nx,:), params.Nx, params.Nx, size(udt,2));
% dTdt = abs(Tdt-Tdt2);
% val = (dTdt>1e-9) & (dTdt>1e-7.*abs(Tdt));
% mean(dTdt(:))
